﻿using System;
using System.Collections.Generic;
using System.Text;
using BarStructures.Deformed;
using SimpleLinearAlgebra;
using MathNet.Numerics.LinearAlgebra;
using BarStructures.Utils;
using MathNet.Numerics;

namespace BarStructures.Analysers
{
    public class FrameAnalyser_XY_CM: Analyser
    {
        private Matrix eigenVectors = null;

        private Vector eigenValues = null;

        public FrameAnalyser_XY_CM(Model model, string name)
            : base(model, name)
        {
            
        }

        public override void Analyse()
        {
            DeformedState ds = new FrameDeformedState();
            ds.AssignToModel(model);

            if (!UseExistingIndexation)
                model.SetAutoSystemIndexes();

            int modelDimension = model.GetNodesCount() * 2;

            //подготовка матриц
            SimpleMatrix stiffness = new SimpleMatrix(modelDimension, modelDimension);
            SimpleMatrix mass = new SimpleMatrix(modelDimension, modelDimension);

            //построение матриц жесткости и масс
            ISimpleMatrix elementMatrix;
            int barsCount, ind1, ind2;
            barsCount = model.GetBarsCount();
            for (int i = 0; i < barsCount; i++)
            {
                Bar bar = model.GetBar(i);
                ind1 = bar.NodeOne.SystemIndex * 2;
                ind2 = bar.NodeTwo.SystemIndex * 2;

                elementMatrix = ElementXYHelper.BarGlobalStiffnessMatrix(bar);
                stiffness.AddSubMatrix(ind1, ind1, elementMatrix.GetSubMatrix(0, 0, 2, 2));
                stiffness.AddSubMatrix(ind1, ind2, elementMatrix.GetSubMatrix(0, 2, 2, 2));
                stiffness.AddSubMatrix(ind2, ind1, elementMatrix.GetSubMatrix(2, 0, 2, 2));
                stiffness.AddSubMatrix(ind2, ind2, elementMatrix.GetSubMatrix(2, 2, 2, 2));

                elementMatrix = ElementXYHelper.BarConsistentMassMatrix(bar);
                mass.AddSubMatrix(ind1, ind1, elementMatrix.GetSubMatrix(0, 0, 2, 2));
                mass.AddSubMatrix(ind1, ind2, elementMatrix.GetSubMatrix(0, 2, 2, 2));
                mass.AddSubMatrix(ind2, ind1, elementMatrix.GetSubMatrix(2, 0, 2, 2));
                mass.AddSubMatrix(ind2, ind2, elementMatrix.GetSubMatrix(2, 2, 2, 2));
            }

            //DbgPrint.SimpleMatrix(mass, "mass_0.txt");
            //DbgPrint.SimpleMatrix(stiffness, "stiffness_0.txt");

            //индексы используемых степеней свободы
            List<int> degreesList = new List<int>();
            int[] degrees;
            for (int i = 0; i < model.GetNodesCount(); i++)
            {
                Node node = model.GetNodeBySystemIndex(i);
                if (!node.Support.Sx)
                    degreesList.Add(i * 2);
                if (!node.Support.Sy)
                    degreesList.Add(i * 2 + 1);
            }
            degrees = degreesList.ToArray();

            ISimpleMatrix K = stiffness.GetSubMatrix(degrees, degrees);
            ISimpleMatrix M = mass.GetSubMatrix(degrees, degrees);

            //DbgPrint.SimpleMatrix(mass, "mass_r.txt");
            //DbgPrint.SimpleMatrix(stiffness, "stiffness_r.txt");

            //решение конденсированной проблемы собственных значений
            Matrix ir_Kcnd = new Matrix(K.GetInternalData());
            Matrix ir_Mcnd = new Matrix(M.GetInternalData());
            Matrix ir = ir_Mcnd.Inverse().Multiply(ir_Kcnd);

            EigenvalueDecomposition evd = new EigenvalueDecomposition(ir);
            eigenVectors = evd.EigenVectors;
            eigenValues = evd.RealEigenvalues;
        }

        public override DeformedState GetDeformedStateForEigenValue(double eigenValue)
        {
            if (eigenValues == null)
                Analyse();

            int i = 0;
            while (i < eigenValues.Length && !Number.AlmostEqual(eigenValue, eigenValues[i]))
                i++;
            if (i < eigenValues.Length)
            {
                DeformedState deformedState = new FrameDeformedState();
                deformedState.AssignToModel(model);
                int counter = 0;
                double dx, dy;
                for (int j = 0; j < model.GetNodesCount(); j++)
                {
                    Node node = model.GetNodeBySystemIndex(j);
                    if (!node.Support.Sx)
                        dx = eigenVectors[counter++, i];
                    else
                        dx = 0;
                    if (!node.Support.Sy)
                        dy = eigenVectors[counter++, i];
                    else
                        dy = 0;
                    NodeDisplacement nodeDisplacement = new NodeDisplacement(node, dx, dy, 0);
                    deformedState.Add(j, nodeDisplacement);
                }

                return deformedState;
            }
            else
                return null;
        }

        public override List<double> GetEigenValuesList()
        {
            if (eigenValues == null)
                Analyse();
            List<double> result = new List<double>();
            for (int i = 0; i < eigenValues.Length; i++)
                result.Add(eigenValues[i]);
            return result;
        }

        public override List<double> GetEigenVecsListForEigenValue(double eigenValue)
        {
            if (eigenValues == null)
                Analyse();

            int i = 0;
            while (i < eigenValues.Length && !Number.AlmostEqual(eigenValue, eigenValues[i]))
                i++;
            List<double> result = new List<double>();
            if (i < eigenValues.Length)
            {
                for (int j = 0; j < eigenVectors.RowCount; j++)
                    result.Add(eigenVectors[j, i]);
            }
            return result;
        }

        public override List<double> GetFullEigenVecsListForEigenValue(double eigenValue)
        {
            throw new NotImplementedException();
        }
    }
}
